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Abstract 

The classical theory of water waves is based on the theory of invis- 
cid flows. However it is important to include viscous effects in some 
applications. Two models are proposed to add dissipative effects in 
the context of the Boussinesq equations, which include the effects of 
weak dispersion and nonlinearity in a shallow water framework. The 
dissipative Boussinesq equations are then integrated numerically. 
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1 Introduction 

Boussinesq equations are widely used in coastal and ocean engineering. 
One example among others is tsunami wave modelling. These equations 
can also be used to model tidal oscillations. Of course, these types of wave 
motion are perfectly described by the Navier-Stokes equations, but currently 
it is impossible to solve fully three-dimensional (3D) models in any significant 
domain. Thus, approximate models such as the Boussinesq equations must 
be used. 

The years 1871 and 1872 were particularly important in the development 
of the Boussinesq equations. It is in 1871 that Valentin Joseph Boussinesq 
received the Poncelet prize from the Academy of Sciences for his work. In 
the Volumes 72 and 73 of the "Comptes Rendus Hebdomadaires des Seances 
de I'Academie des Sciences" , which cover respectively the six-month periods 
January- June 1871 and July-December 1871, there are several contributions 
of Boussinesq. On June 19, 1871, Boussinesq presents the now famous note 
on the solitary wave entitled "Theorie de I'intumescence liquide appelee onde 
solitaire ou de translation, se propageant dans un canal rectangulaire" (72, 
pp. 755-759), which will be extended later in the note entitled "Theorie 
generale des mouvements qui sont propages dans un canal rectangulaire hor- 
izontal" (73, pp. 256-260). Saint- Venant presents a couple of notes of Boussi- 
nesq entitled "Sur le mouvement permanent varie de I'eau dans les tuyaux 
de conduite et dans les canaux decouverts" (73, pp. 34-38 and pp. 101- 
105). Saint- Venant himself publishes a couple of notes entitled "Theorie du 
mouvement non permanent des eaux, avec application aux crues des rivieres 
et a I'introduction des marees dans leur lit" (73, pp. 147-154 and pp. 237- 
240). All these notes deal with shallow-water theory. On November 13, 1871, 
Boussinesq submits a paper entitled "Theorie des ondes et des remous qui se 
propagent le long d'un canal rectangulaire horizontal, en communiquant au 
liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface 
au fond", which will be published in 1872 in the Journal de Mathematiques 
Pures et Appliquees (17, pp. 55-108). 

Boussinesq [1872, 1871] included dispersive effects for the first time in the 
Saint- Venant equations [de Saint- Venant, 1871]. One should mention that 
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Boussinesq's derivation was restricted to 1 + 1 dimensions (x and t) and a hor- 
izontal bottom. Boussinesq equations contain more physics than the Saint- 
Venant equations but at the same time they are more comphcated from the 
mathematical and numerical point of views. These equations possess a hyper- 
bolic structure (the same as in the nonlinear shallow-water equations) com- 
bined with high-order derivatives to model wave dispersion. There have been 
a lot of further developments of these equations like in Madsen and Schaffer 
[1998], Nwogu [1993], Peregrine [1967], Wei et al. [1995]. 

Let us outline the physical assumptions. The Boussinesq equations are in- 
tended to describe the irrotational motion of an incompressible homogeneous 
inviscid fluid in the long wave limit. The goal of this type of modelling is to 
reduce 3D problems to two-dimensional (2D) ones. This is done by assum- 
ing a polynomial (usually linear) vertical distribution of the flow field, while 
taking into account non-hydrostatic effects. This is the principal physical 
difference with the nonlinear shallow-water (NSW) equations. 

There are a lot of forms of the Boussinesq equations. This diversity is due 
to different possibilities in the choice of the velocity variable. In most cases 
one chooses the velocity at an arbitrary water level or the depth-averaged 
velocity vector. The resulting model performance is highly sensitive to linear 
dispersion properties. The right choice of the velocity variable can signifi- 
cantly improve the propagation of moderately long waves. A good review 
is given by Kirby [2003]. There is another technique used by Bona et al. 
[2002]. Formally, one can transform higher-order terms by invoking lower- 
order asymptotic relations. It provides an elegant way to improve the proper- 
ties of the linear dispersion relation and it gives a quite general mathematical 
framework to study these systems. 

The main purpose of this article is to include dissipative effects in the 
Boussinesq equations. It is well-known that the effect of viscosity on free 
oscillatory waves on deep water was studied by Lamb [1932]. What is less 
known is that Boussinesq himself studied this effect as well. Boussinesq 
wrote three related papers in 1895 in the "Comptes Rendus Hebdomadaires 
des Seances de I'Academie des Sciences": (i) "Sur I'extinction graduelle de 
la houle de mer aux grandes distances de son lieu de production : formation 
des equations du probleme" (120, pp. 1381-1386), (ii) "Lois de I'extinction 
de la houle en haute mer" (121, pp. 15-20), (iii) "Sur la maniere dont se 
regularise au loin, en s'y reduisant a une houle simple, toute agitation confuse 
mais periodique des fiots" (121, pp. 85-88). It should be pointed out that the 
famous treatise on hydrodynamics by Lamb has six editions. The paragraphs 
on wave damping are not present in the first edition (1879) while they are 
present in the third edition (1906). The authors did not have access to the 
second edition (1895), so it is possible that Boussinesq and Lamb published 



1 Introduction 



4 



similar results at the same time. Indeed Lamb derived the decay rate of the 
linear wave amplitude in two different ways: in paragraph 348 of the sixth 
edition by a dissipation calculation (this is also what Boussinesq [1895] did) 
and in paragraph 349 by a direct calculation based on the linearized Navier- 
Stokes equations. Let a denote the wave amplitude, v the kinematic viscosity 
of the fluid and k the wavenumber of the decaying wave. Boussinesq (see Eq. 
(12) in Boussinesq [1895]) and Lamb both showed that 



Equation (1) leads to the classical law for viscous decay of waves of amplitude 
a, namely a ~ exp(— 2z//i;^t) (see Eq. (13) in Boussinesq [1895] after a few 
calculations). 

In the present paper, we use two different models for dissipation and de- 
rive corresponding systems of long-wave equations. There are several meth- 
ods to derive the Boussinesq equations but the resulting equations are not 
the same. So one expects the solutions to be different. We will investigate 
numerically whether corresponding solutions remain close or not. 

One may ask why dissipation is needed in Boussinesq equations. First 
of all, real world liquids are viscous. This physical effect is "translated" in 
the language of partial differential equations by dissipative terms (e.g. the 
Laplacian in the Navier-Stokes equations). So, it is natural to have analogous 
terms in the long wave limit. In other words, a non-dissipative model means 
that there is no energy loss, which is not pertinent from a physical point of 
view, since any flow is accompanied by energy dissipation. 

Let us mention an earlier numerical and experimental study by Bona et al. 
[1981]. They pointed out the importance of dissipative effects for accurate 
long wave modelling. In the "Resume" section one can read 

[...] it was found that the inclusion of a dissipative term was 
much more important than the inclusion of the nonlinear term, 
although the inclusion of the nonlinear term was undoubtedly 
beneficial in describing the observations [...]. 

The complexity of the mathematical equations due to the inclusion of this 
term is negligible compared to the benefit of a better physical description. 

Let us consider the incompressible Navier-Stokes (N-S) equations for a 
Newtonian fluid: 





V ■ u* = 



(9u* 



+ u* ■ Vu* 



+ z/Au* + — , 
P 



dt* 



P 
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where u*{x,y,z,t) = {u* ,v* ,w*){x* ,y* , z* ,t*) is the fluid velocity vector, p* 
the pressure, F* the body force vector, p the constant fluid density and z/ the 
kinematic viscosity. 

Switching to dimensionless variables by introducing a characteristic veloc- 
ity U, a characteristic length L and a characteristic pressure pU"^, neglecting 
body forces^ in this discussion, the N-S equations become 

V-u = 0, 

9u ^ „ 1 ^ 
— + u ■ Vu = -Vp + ^Au, 

ot Re 

where Re is the well-known dimensionless parameter known as the Reynolds 
number and deflned as 

Finertia U L 



Re 



F ■ V 

^ VISCOUS 



From a physical point of view the Reynolds number is a measure of the 
relative importance of inertial forces compared to viscous effects. For typ- 
ical tsunami propagation applications the characteristic particle velocity U 
is about 5 cm/s and the characteristic wave amplitude, which we use here 
as characteristic length scale, is about 1 m. The kinematic viscosity u de- 
pends on the temperature but its order of magnitude for water is 10~^ m^/s. 
Considering that as the tsunami approaches the coast both the particle ve- 
locity and the wave amplitude increase, one can write that the corresponding 
Reynolds number is of the order of 10^ or 10^. This simple estimate clearly 
shows that the flow is turbulent (as many other flows in nature). 

It is a common practice in fluid dynamics (addition of an "eddy viscosity" 
into the governing equations for Large Eddy Simulations^) to ignore the 
small-scale vortices when one is only interested in large-scale motion. It can 
signiflcantly simplify computational and modelling aspects. So the inclusion 
of dissipation can be viewed as the simplest way to take into account the 
turbulence. 

There are several authors [Dias et al., 2007, Longuet-Higgins, 1992, Ruvinsky et al., 
1991, Skandrani et al., 1996, Spivak et al., 2002, Tuck, 1974] who included 
dissipation due to viscosity in potential flow solutions and there are also 
authors [Heitner and Housner, 1970, Kennedy et al., 2000, Zelt, 1991] who 
already included in Boussinesq models ad-hoc dissipative terms into momen- 
tum conservation equations in order to model wave breaking. Modelling this 
effect is not the primary goal of the present paper, since the flow is no longer 



^The presence or absence of body forces is not important for discussing viscous effects. 
^Boussinesq himself introduced the concept of eddy viscosity in his famous 680 page 
paper entitled "Essai sur la thcorie des eaux courantes" [Boussinesq, 1877]. 
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irrotational after wave breaking. Strictly speaking tlie Boussinesq equations 
can no longer be valid at this stage. Nevertheless scientists and engineers 
continue to use these equations even to model the run-up on the beach. In 
our approach a suitable choice of the eddy viscosity, which is a function of 
both space and time, can model wave breaking at least as well as in the 
articles cited above. 



2 Derivation of the Boussinesq equations 




Figure 1: Sketch of the fluid domain 

In order to derive the Boussinesq equations, we begin with the full water- 
wave problem. A Cartesian coordinate system (x', y', z') is used, with the x'— 
and y'—axis along the still water level and the z'— axis pointing vertically 
upwards. Let fit be the fluid domain in M.^ which is occupied by an inviscid 
and incompressible fluid. The subscript t underlines the fact that the domain 
varies with time and is not known a priori. The domain Qf is bounded 
below by the seabed z' = —h'{x',y',t') and above by the free surface z' = 
rj' {x' , y' , t') . In this section we choose the domain Qt to be unbounded in 
the horizontal directions in order to avoid the discussion on lateral boundary 
conditions. The reason is twofold. First of all, the choice of the boundary 
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value problem (BVP) (e.g. generating and/or absorbing boundary, wall, run- 
up on a beach) depends on the application under consideration and secondly, 
the question of the well-posedness of the BVP for the Boussinesq equations is 
essentially open. Primes stand for dimensional variables. A typical sketch of 
the domain Qt is given in Figure 1. If the flow is assumed to be irrotational 
one can introduce the velocity potential (p' defined by 

' \dx' ' dy'' dz' J ' 

where u' denotes the velocity field. Then we write down the following system 
of equations for potential fiow theory in the presence of a free surface: 

A'0' = 0, (x', y', z') G fit = X l-h', r]'], 



4, + l\V4>f + gr]' = 0, z' = rj', (2) 

+ h[, + V>' ■ Vh' = 0, z' = -h\ 

where g denotes the acceleration due to gravity (surface tension effects are 
usually neglected for long- wave applications). It has been assumed implicitly 
that the free surface is a graph and that the pressure is constant on the free 
surface (no forcing). Moreover we assume that the total water depth remains 
positive, i.e. 77' + /i' > (there is no dry zone). 

As written, this system of equations does not contain any dissipation. 
Thus, we complete the free-surface dynamic boundary condition (2) by adding 
a dissipative term to account for the viscous effects^: 

In this work we investigate two models for the dissipative term D'^,. For 
simplicity, one can choose a constant dissipation model (referred hereafter as 
Model I) which is often used (e.g. [Jiang et al., 1996]): 

Model I: := (3) 

■^Dias et al. [2007], who considered deep-water waves, pointed out that a viscous cor- 
rection should also be added to the kinematic boundary condition if one takes into account 
the vortical component of the velocity. This correction was recently added in finite depth 
as well [Dutykh and Dias, 2007a]. A boundary-layer correction at the bottom was also 
included. 
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There is a more physically realistic dissipation model which is obtained 
upon balancing of normal stress on the free surface (e.g. [Dias et al., 2007, 
Ruvinsky et al., 1991, Zhang and Vinals, 1997]): 

Model II: := (520',,^,. (4) 

The derivation of Boussinesq equations is more transparent when one 
works with scaled variables. Let us introduce the following independent and 
dependent non-dimensional variables: 



ho ao gaoi 

where ho, i and ao denote a characteristic water depth, wavelength and wave 
amplitude, respectively. 

After this change of variables, the set of equations becomes 

/i^(0xx + 0ot) + = 0' {x,y,z)ent, (5) 



-'z 



/i^?7t + e/i^V0 ■ V?7, z = er], (6) 



fi% + ^efi^\V<j)f + ^e<Pl + fi^r] + eD^ = 0, z = er] (7) 

<Pz + + /iV0 ■ V/i = 0, z = -h, (8) 

where e and fi are the classical nonlinearity and frequency dispersion param- 
eters defined by 

ao ho 
ho i 

In these equations and hereafter the symbol V denotes the horizontal gradi- 
ent: 

d d^^ 



The dissipative term is given by the chosen model (3) or (4): 



Model I: = ^0, Model II: = ^ 
Ri R2 



■'ZZ1 



2 Derivation of the Boussinesq equations 



9 



parameter value 

Acceleration due to gravity g, m/s^ 10 

Amplitude oq, m 1 

Wave length i, km 100 

Water depth h^, km 4 

Kinematic viscosity 6, m^/s 10~^ 



Table 1: Typical values of characteristic parameters in tsunami applications 



where the following dimensionless numbers have been introduced: 

^ 1 / gapi \ ^ _ 1 / gapi \ 

^' ^i\hly/gho)' ^' ^2\Vgho) 

From this dimensional analysis, one can conclude that the dimension of the 
coefficient 6i is [s~^] and that of §2 is [m^s~^]. Thus, it is natural to call the 
first coefficient viscous frequency (since it has the dimensions of a frequency) 
and the second one kinematic viscosity (by analogy with the N-S equations). 

It is interesting to estimate R2, since we know how to relate the value of 
62 to the kinematic viscosity of water. Typical parameters which are used 
in tsunami wave modelling are given in Table 1. For these parameters R2 = 
5 X 10^ and fi^ = 1.6 x 10"'^. The ratio between inertial forces and viscous 
forces is ^efi'^\V(f)\'^ /e\D^\. Its order of magnitude is fi^R2, that is 8 x 10®. It 
clearly shows that the flow is turbulent and eddy-viscosity type approaches 
should be used. It means that, at zeroth-order approximation, the main 
effect of turbulence is energy dissipation. Thus, one needs to increase the 
importance of viscous terms in the governing equations in order to account 
for turbulent dissipation. 

As an example, we refer one more time to the work by Bona et al. [1981]. 
They modeled long wave propagation by using a modifled dissipative Korteweg- 
de Vries equation: 

3 1 

Vt + Vx + TiVVx - l^rjxx - -Vxxt = 0. (9) 
2 

In numerical computations the authors took the coefficient fi = 0.014. This 
value gave good agreement with laboratory data. 

From now on, we will use the notation Ui := 1/Ri. This will allow us to 
unify the physical origin of the numbers Ri with the eddy- viscosity approach. 
In other words, for the sake of convenience, we will "forget" about the origin 
of these coefficients, because their values can be given by other physical 
considerations. 
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2.1 Asymptotic expansion 

Consider a formal asymptotic expansion of tlie velocity potential in 
powers of the small parameter /z^: 

(t) = (t>0 + /iVl + /iV2 + • • ■ ■ (10) 

Then substitute this expansion into the continuity equation (5) and the 
boundary conditions. After substitution, the Laplace equation becomes 

/x2(V2(/.o + + /i'V'</.2 + ...)+ <Pozz + + /xV2.. + . . . = 0. 

Collecting the same order terms yields the following equations in the domain 

/ : <Pozz = 0, (11) 
/i': 0i.. + VVo = O, (12) 
02.. + VVi = O. (13) 

Performing the same computation for the bottom boundary condition yields 
the following relations aX z = —h: 

: <^o. = 0, (14) 

/x': (/.I, + i/ii + V0O ■ V/i = 0, (15) 
/i^ : 02. + V01 ■ V/i = 0. (16) 

From equation (11) and the boundary condition (14) one immediately con- 
cludes that 

Let us define the horizontal velocity vector 

u(x,y,t) := V0O, u=(u,t;)^. 

The expansion of Laplace equation in powers of /i^ gives recurrence rela- 
tions between 0o, 0i, 02, etc. Using (12) one can express 0i in terms of the 
derivatives of 0o: 

012. = -V ■ u. 
Integrating once with respect to z yields 

01. = -zV ■ u + Ci(x, y, t). 

The unknown function Ci{x,y,t) can be found by using condition (15): 

0iz = —{z + h)V ■ u ht — u ■ V/i, 
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and integrating one more time witli respect to z gives the expression for 0i: 
= -^(z + /i)V-u-^ (-/if + u- V/z ) . (17) 



Now we will determine 02- For this purpose we use equation (13): 
"y^zz = \{z + /i)'V'(V ■ u) + [{h + z)V'/i + |V/i|') V ■ u+ 



2{h + 2)V/i ■ V(V ■ u) + 2 ( ^V^/it + V^(u ■ V/i) 



Integrating twice with respect to z and using the bottom boundary con- 
dition (16) yields the following expression for 

02 = ^ih + z)^V'(V ■ u) + {-{z + /i)3v'/i + \z^ |V/irt V ■ u 

+ ^(z + Kf^h ■ V(V ■ u) + ^(-V'/it + V'(u ■ V/i)) 

-2/i(^V^(-/it + u- V/i) + V/i- V(-/it + u- V/i) - |V/i|^V-u). (19) 

Remark: In these equations one finds the term {l/e)ht due to the moving 
bathymetry. We would like to emphasize that this term is 0(1), since in 
problems of wave generation by a moving bottom the bathymetry h{x, y, t) 
has the following special form in dimensionless variables: 

h{x, y, t) := ho{x, y) - eC{x, y, t), (20) 

where ho{x,y) is the static seabed and ({x,y,t) is the dynamic component 
due to a seismic event or a landslide (see for example Dutykh and Dias 
[2007b] for a practical algorithm constructing ({x,y,t) in the absence of 
a dynamic source model). The amplitude of the bottom motion has to be 
of the same order of magnitude as the resulting waves, since we assume the 
fluid to be inviscid and incompressible. Thus {l/e)ht = —Q = 0(1). 

In the present study we restrict our attention to dispersion terms up to 
order 0{fi'^). We will also assume that the Ursell- Stokes number [Ursell, 
1953] is 0(1): 

5:= 4 = 0(1). 

This assumption implies that terms of order 0(£^) and 0(e/i^) must be ne- 
glected, since 
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Of course, it is possible to obtain high-order Boussinesq equations. We de- 
cided not to take this research direction. For high-order asymptotic expan- 
sions we refer to Madsen and Schaffer [1998], Wei et al. [1995]. Recently, 
[Benoit, 2006] performed a comparative study between fully- nonhnear equa- 
tions [Wei et al., 1995] and Boussinesq equations with optimized dispersion 
relation [Nwogu, 1993]. No substantial difference was revealed. 

Now, we are ready to derive dissipative Boussinesq equations in their 
simplest form. First of all, we substitute the asymptotic expansion (10) into 
the kinematic free-surface boundary condition (6): 



+ /i^^i, + /02. = fi^Vt + eyU^V^o ■V7] + 0{e^ + + yu^), z = er]. (21) 



The first term on the left hand side is equal to zero because of Eq. (14). 

Using expressions (17) and (19) one can evaluate (piz and (j)2z on the free 
surface: 



-{h + er])V ■u--ht-u-Vh 



=eri 



•J2z 



■ + ■ ■ u) + /i ^ + 1 v/irj V ■ u 

1 1 
- ~ h-^ht ■ Vh + 0{e). 

Substituting these expressions into (21) and retaining only terms of order 
0{e + fi"^) yields the free-surface elevation equation: 

Vt + ^-(ih + eri)u) = -( 1 + ^h^V^ + fx^hVh ■ -ht + ii'^V^V ■ n) 

\ 2 / e 6 

+ /i^/i (^Vh ■ V(V ■ u) + Qv^/i + \Vhf^ V ■ . 

The equation for the evolution of the velocity field is derived similarly 
from the dynamic boundary condition (7). This derivation will depend on 
the selected dissipation model. For both models one has to evaluate 0i, (pu 
and (pizz along the free surface z = erj and then substitute the expressions 
into the asymptotic form of (7): 

AiVoi + /i^0it + |V0o|^ + + £i^2fJ''^(pizz = 0{e^ + efi'^ + /i^), 

where, as an example, dissipative terms are given according to the second 
model. After performing all these operations one can write down the follow- 
ing equations: 

Model I: (Pot + fu^ + r/ + z/i^0o - ^/^'V ■ u - i^h^V ■ = 0, 
Model II: (pot + fu^ + - z/s^V ■ u - ^/i^V • u< = 0. 
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The last step consists in differentiating the above equations with respect to 
the horizontal coordinates in order to obtain equations for the evolution of 
the velocity. We also perform some minor transformations using the fact that 
the vector u is a gradient by definition, so we have the obvious relation 

du dv 
dy dx 

The resulting Boussinesq equations for the first and second dissipation 
models, respectively, are given below: 

r]t + V-({h + e7])n) = - (l + y ^'V' + fi^hWh ■ \ + /i'^ V'(V ■ u) 



(22) 



+fi^h (^hVh ■ V(V ■ u) + Qv^/i + \Vhf'^ V ■ u 

Model I: + ^eVu^ + V?7 + uiSu = ^euiV^h'^V ■ u) + \^'^V{h'^V ■ u^y 
Model II: VLt + \eWvL^ + Vr/ = ez/aV^u + \fi^V{h^W ■ ut). (24) 

3 Analysis of the linear dispersion relations 

For simplicity, we will consider in this section only 2D problems. The 
generalization to 3D problems is straightforward and does not change the 
analysis. 

3.1 Linearization of the full potential flow equations 
with dissipation 

First we write down the linearization of the full potential flow equations 
in dimensional form, after dropping the primes: 

A0 = O, (x,^) G M X [-/i,0], (25) 

0. = Vu z = 0, (26) 

0i + + Z^^ = 0, z = 0, (27) 

0, = 0, z = -h. (28) 

Remark: In this section the water layer is assumed to be of uniform 
depth, so h = const. 
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As above the term D^p depends on the selected dissipation model and is 
equal to 5i0 or 62(pzz- The next step consists in choosing a special form of 
solutions: 

z, t) = ifoe'^''''-'^'^if{z), 7]{x, t) = r/oe^('=^-"*\ (29) 

where (po and rjQ are constants. Substituting this form of solutions into equa- 
tions (25), (26) and (28) yields the following boundary value problem for an 
ordinary differential equation: 

<^"(z) - fcV(^) = 0, ze[-h,ol 

<^'(0) = ^i-iuj), ^\-h) = 0. 
Straightforward computations give the solution to this problem: 

The dispersion relation can be thought as a necessary condition for so- 
lutions of the form (29) to exist. The problem is that uj and k cannot be 
arbitrary. We obtain the required relation uj = uj{k), which is called the 
dispersion relation, after substituting this solution into (27). 

When the dissipative term is chosen according to model I (3), = 6i(p 
and the dispersion relation is given implicitly by 

tu^ -|- i5iuj — gktanh{kh) = 0, 

or in explicit form by 



Lj = ±\l gktanhikh) - ^ (30) 

For the second dissipation model (4) one obtains the following relation: 
cj^ -|- i62Ujk'^ — gkta.nh.{kh) = 0. 
One can easily solve this quadratic equation for as a function of k: 



UJ = ±^gktanh{kh) - - Y^^- (31) 

If 6i,2 = one easily recognizes the dispersion relation of the classical 
water-wave problem: 

u = ±^/gktanh{kh). (32) 
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Remark: It is important to have the property lmuj{k) < 0,V/c in order to 
avoid the exponential growth of certain wavelengths, since 



For our analysis it is more interesting to look at the phase speed which 
is defined as 




The phase velocity is directly connected to the speed of wave propagation and 
is extremely important for accurate tsunami modelling since tsunami arrival 
time obviously depends on the propagation speed. The expressions for the 
phase velocity are obtained from the corresponding dispersion relations (30) 
and (31): 



It can be shown that in order to keep the phase velocity unchanged by the 
addition of dissipation, similar dissipative terms must be included in both 
the kinematic and the dynamic boundary conditions [Dias et al., 2007]. 

3.2 Dissipative Boussinesq equations 

The analysis of the dispersion relation is even more straightforward for 
Boussinesq equations. In order to be coherent with the previous subsection, 
we switch to dimensional variables. As usual we begin with the (1 + 1)D 
linearized equations: 



,i(kx—uj{k)t) liaLu(k)t i(kx—'ReLj{k)t) 




(33) 




(34) 



r]t + hu^ = —u, 



'XXX J 



Model II 



Model I 



Ut + gVx = 52Uxx + \h^Uxxt- 



Now we substitute a special ansatz in these equations: 



,i(kx—u!t) 



U = UqC' 



^i{kx—Lot) 
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where rjo and uq are constants. In the case of the first model, one obtains 
the following homogeneous system of linear equations: 

{-iuj)r]o + ikh ^1 + ^(kh)'^^ uq =0, 
gikrjo + (^-iu + 5i + jikhf - y (^/^)^^ ^^o = 0. 

This system admits nontrivial solutions if its determinant is equal to zero. It 
gives the required dispersion relation: 

A similar relation is found for the second model: 



The corresponding phase velocities are given by 



^^M\^-(^^% (35) 



pb y^"\i + Ukh)y \2k) 2k 



(2) _ / , / i + _ ( S^k V tS^k 

""pb -\l9^[i+ukhy) y2 + {khy) 2 + {khy ^^^^ 



3.3 Discussion 

Let us now provide a discussion on the dispersion relations. The real and 
imaginary parts of the phase velocities (33)-(36) for the full and long wave 
linearized equations are shown graphically on Figures 2-7. In this example 
the parameters are given by Si = 0.14,^2 = 0.14. Together with the dissi- 
pative models we also plotted for comparison the well-known phase velocity 
corresponding to the full conservative (linearized) water-wave problem: 



, , , tanh(kh) 

First of all, one can see that dissipation is very selective, as is often the 
case in physics. Clearly, the first dissipation model prefers very long waves, 
while the second model dissipates essentially short waves. Moreover one can 
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Figure 2: Dissipation model I. Real part of the phase velocity. 
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Figure 3: Dissipation model I. Same as Figure 2 with a zoom on long waves. 



3.3 Discussion 



18 



Dispersion relation. Imaginary part 



-0.02 



-0.04 



-0.06 



-0.08 



-0.1 



-0.12 



-0.14 



0.05 



s 






\ 

\ 


Full equations with dissipation-1 . Branch - 

Full equations with dissipation-1 . Branch + 




\ 

\ 

\ 


Dissipative Boussinesq-1 












\ 
1 






1 
1 






1 
1 






1 
1 

. . . . I 






/ 
/ 

/ 






/ 

/ 

/ 







0.1 
kh 



0.15 



0.2 



Figure 4: Dissipation model I. Imaginary part of the frequency. 
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Figure 5: Dissipation model II. Real part of the phase velocity. 
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Figure 6: Dissipation model II. Imaginary part of the frequency. 
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Figure 7: Dissipation model II. Same as Figure 6 with a zoom on long waves. 
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see from the expressions (33), (35) that the phase velocity has a l//c singular 
behaviour in the vicinity of kh = (in the long wave limit). Furthermore, it 
can be clearly seen in Figure 3 that very long linear waves are not advected 
in the first dissipation model, since the real part of their phase velocity is 
identically equal to zero. 

That is why we suggest to make use of the second model in applications 
involving very long waves such as tsunamis. 

On the other hand we would like to point out that the second model 
admits a critical wavenumber kc such that the phase velocity (34) becomes 
purely imaginary with negative imaginary part. From a physical point of 
view it means that the waves shorter than kc are not advected, but only dis- 
sipated. When one switches to the Boussinesq approximation, this property 
disappears for physically realistic values of the parameters g, h and 62 (see 



Let us clarify this situation. The qualitative behaviour of the phase veloc- 
ity c^pjj (see equation (36)) depends on the roots of the following polynomial 
equation: 



This equation does not have real roots since ^ -C 1. 

4 Alternative version of the Boussinesq equa- 
tions 

In this section we give an alternative derivation of Boussinesq equa- 
tions. We use another classical method for deriving Boussinesq-type equa- 
tions [Benjamin, 1974, Peregrine, 1972, Whitham, 1999], which provides 
slightly different governing equations. Namely, the hyperbolic structure is 
the same, but the dispersive terms differ. In numerical simulations we sug- 
gest to use this system of equations. 

The derivation follows closely the paper by Madsen and Schaffer [1998]. 
The main differences are that we neglect the terms of order 0(yU,^), take in 
account a moving bathymetry and, of course, dissipative effects which are 
modelled this time according to model II (4) because, in our opinion, this 
model is more appropriate for long wave applications. Anyhow, the derivation 
process can be performed in a similar fashion for model I (3). 



Table 1). 




{khY + 12 = 0. 
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4.1 Derivation of the equations 

The starting point is the same: equations (5), (6), (7) and (8). This time 
the procedure begins with representing the velocity potential y, z, t) as 
a formal expansion in powers of z rather than of /x^: 

oo 

y, z,t) = ^ z'^(pn{x, y, t). (37) 

n=0 

We would like to emphasize that this expansion is only formal and no con- 
vergence result is needed. In other words, it is just convenient to use this 
notation in asymptotic expansions but in practice, seldom more than four 
terms are used. It is not necessary to justify the convergence of the sum with 
three or four terms. 

When we substitute the expansion (37) into Laplace equation (5), we 
have an infinite polynomial in z. Requiring that formally satisfies Laplace 
equation implies that the coefficients of each power of z vanish (since the 
right-hand side is identically zero). This leads to the classical recurrence 
relation 

0n,+2(a;, t) = --, — — TTT — rV^0„(x, t), n = 0, 1, 2, . . . 
(n + l)(n + 2) 

Using this relation one can eliminate all but two unknown functions in (37): 

oo /In 2n+l 

The following notation is introduced: 

uo := u(a;,?/,0,t), Wq := \w{x,y,^,t). 



It is straightforward to find the relations between Uq, Wq and 0o; 0i if one 

d_ 

dz 



remembers that (u, w) = (V, ■§-)(t>'- 



Uo = V0O, Wo = ^01. 

Using the definition of the velocity potential cj) one can express the velocity 
field in terms of uq, Wq. 

oo / 2n 2n+l \ 
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■JO / 2n+l In \ 

- - E(-l)" (- (^"^"^^V-lV . u„) + ^.-"V-..„) . 

These formulas are exact but not practical. In the present work we neglect 
the terms of order 0(/i^) and higher. In this asymptotic framework the above 
formulas become much simpler: 

= 00 + ^01 - ^ (V^^o + |v20i) + 0(/i^), (38) 

2 2 

u = uo + z/i'Vi/;o-^V(V-uo) + 0(/i^), (39) 

w = i?wq- zi?^ ■\iq\0{ii^). (40) 

In order to establish the relation between and uq one uses the bot- 
tom kinematic boundary condition (8), which has the following form after 
substituting the asymptotic expansions (38), (39), (40) in it: 



ht + e^h ■ (^uo - hi^'^Vwo - V(V • uq)) 

2V72^V7 N ^^^^ 



+ e{wo- y/iV2(V ■ uo) - ^V^t^o) + 0(/) = 0. (41) 

In order to obtain the expression of Wq in terms of uq one introduces one 
more expansion: 

Wo{x,y,t) = wf\x,y,t) + ij,'^w^^\x,y,t) + .... (42) 

We insert this expansion into the asymptotic bottom boundary condition 



(41). This leads to the following explicit expressions for w^l^^ and Wq^^' 



4'^ = y(v/^-v(v-uo) + ^v2(v-uo; 



- h (vh ■ ^Vht + Vh ■ V(V ■ (huo)) + |(^V'/it + V'(V ■ (/iuq)))^ . 

Substituting these expansions into (42) and performing some simplifications 
yields the required relation between uq and Wq: 



Wq = -^ht - V ■ (huo) 



/iV ■ + y V(V ■ (hu,)) - -V(V ■ Uo)) + 0(/). (43) 
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Now one can eliminate the vertical velocity wo since one has its expression 
(43) in terms of uq. Equations (39)-(40) become 

2 2 

u = uo - z^Vht - i2'(^zV{V ■ {huo)) + y V(V ■ uo)) + 0(/), (44) 

2 

w = -^ht - /i^ (V ■ (huo) + zV ■ Uo) + 0(/i^). (45) 

In this work we apply a trick due to Nwogu [1993] . Namely, we introduce 
a new velocity variable u^ defined at an arbitrary water level Za = —ah. 
Technically this change of variables is done as follows. First we evaluate (44) 
ai z = Za, which gives the connection between uq and Uq,: 

2 2 

u„ = Uo - Za^Vht - /i'(^«V(V ■ (huo)) + ^V(V ■ uo)) + 0(/). 

Using the standard technics of inversion one can rewrite the last expression 
as an asymptotic formula for Uq in terms of u^: 

2 2 

Uo = u, + za^ht + fj.^ (^«V(V • {hua)) + ^V(V ■ u,)) + 0(/). (46) 

Remark: Behind this change of variables there is one subtlety which is 
generally hushed up in the literature. In fact, the wave motion is assumed to 
be irrotational since we use the potential fiow formulation (5), (6), (7), (8) of 
the water-wave problem. By construction rot(u, w) = when u and w are 
computed according to (44), (45) or, in other words, in terms of the variable 
Uq. When one turns to the velocity variable defined at an arbitrary level, 
one can improve the linear dispersion relation and this is important for wave 
modelling. But on the other hand, one loses the property that the fiow is 
irrotational. That is to say, a direct computation shows that rot(u, w) ^ 
when u and w are expressed in terms of the variable u^. The purpose of 
this remark is simply to inform the reader about the price to be paid while 
improving the dispersion relation properties. It seems that this point is not 
clearly mentioned in the literature on this topic. 

Let us now derive the Boussinesq equations. There are two different 
methods to obtain the free-surface elevation equation. The first method 
consists in integrating the continuity equation (5) over the depth and then 
use the kinematic free-surface and bottom boundary conditions. The second 
way is more straightforward. It consists in using directly the kinematic free- 
surface boundary condition (6): 



r]t + 5V0 ■ Vt] - —(j), = 0. 
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Then one can substitute (38) into (6) and perform several simplifications. 
Neglecting all terms of order 0{e^ + + /x^) yields the following equation^ 

2 73 

r/i + V ■ ((/i + £r/)uo) + ■ (/iV(V ■ (/iuq)) - y V(V ■ uq)) = 

= + yV-(Wc*). 

Recall that ({x,y,t) is defined according to (20). When the bathymetry is 
static, C = 0. We prefer to introduce this function in order to eliminate the 
division by e in the source terms since this division can give the impression 
that stiff source terms are present in our problem, which is not the case. 

In order to be able to optimize the dispersion relation properties, we 
switch to the variable u^. Technically it is done by using the relation (46) 
between uq and Uq,. The result is given below: 

/ h h h'^ \ 

r/t+V-((/i+er/)u«)+/i'V-(^/i(2„ + -)V(V-(/iu„)) + -(2;2-_)v(V-u,)j = 

= Ci + yu'V- (/i(^, + |)V0). (47) 

As above, the equation for the horizontal velocity field is derived from the 
dynamic free-surface boundary condition (7). It is done exactly as in section 
2 and we do not insist on this point: 

uot + |V |uo|' + Vr/ - e6V^uo = 0. 

Switching to the variable Uq, yields the following governing equation: 

2 

Uat + |v|u,|2 + Vr/ + /i2(2;«V(V ■ (hn^)) + y V(V ■ u,))^ = 

= e6Auo, + fi\zaVCt)t. (48) 

In several numerical methods it can be advantageous to rewrite the system 
(47), (48) in vector form: 

Vt + /i'L(U)t + V ■ F(U) + fi^V ■ P(U) = S(x, y, t) + eSV ■ (DVU)), 

where 

('^\ ^ X. dF dG 
U:= UJ , V-F := — + —, 

^We already discussed this point on page 11. In this section we also assume that the 
Stokes- Ursell number S is of order 0(1). 
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4.2 Improvement of the linear dispersion relations 

As said above, the idea of using one free parameter a G [0, 1] to optimize 
the linear dispersion relation properties appears to have been proposed first 
by Nwogu [1993]. 

The idea of manipulating the dispersion relation was well-known before 
1993. See for example Madsen et al. [1991], Murray [1989]. But these authors 
started with a desired dispersion relation and artificially added extra terms 
to the momentum equation in order to produce the desired characteristics. 
We prefer to follow the ideas of Nwogu [1993]. 

Remark: When one plays with the dispersion relation it is important to 
remember that the resulting problem must be well-posed, at least linearly. 
We refer to Bona et al. [2002] as a general reference on this topic. Usually 
Boussinesq-type models with good dispersion characteristics are linearly well- 
posed as well. 

In order to look for an optimal value of a we will drop dissipative terms. 
Indeed we want to concentrate our attention on the propagation properties 
which are more important. 

The choice for the parameter a depends on the optimization criterion. 
In the present work we choose a by comparing the coefficients in the Taylor 
expansions of the phase velocity in the vicinity of kh = 0, which corresponds 
to the long- wave limit. Another possibility is to match the dispersion relation 
of the full linearized equations (32) in the least square sense. One can also 
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use Pade approximants [Witting, 1984] since rational functions have better 
approximation properties than polynomials. 

We briefly describe the procedure. First of all one has to obtain the phase 
velocity of the linearized, no n- viscous, Boussinesq equations (47)- (48). The 
result is 

gh i_a(|_i)(fc/i)2 3^ ^ 6 ^ ' ' ^ 

(49) 

On the other hand one can write down the phase velocity of the full linearized 
equations (32): 

If one insists on the dispersion relation (49) to be exact up to order O {{khY) 
one immediately obtains an equation for aopt- 



^opt) ^ 

6 =15^"°^* = ^- — 



0.55. 



We suggest using this value of a in numerical computations. 

4.3 Bottom friction 

In this subsection, one switches back to dimensional variables. It is a 
common practice in hydraulics engineering to take into account the effect of 
bottom friction or bottom rugosity. In the Boussinesq and nonlinear shallow 
water equations there is also a possibility to include some kind of empirical 
terms to model these physical effects. From the mathematical and especially 
numerical viewpoints these terms do not add any complexity, since they have 
the form of source terms that do not involve differential operators. So it is 
highly recommended to introduce these source terms in numerical models. 

There is no unique bottom friction law. Most frequently, Chezy and 
Darcy-Weisbach laws are used. Both laws have similar structures. We give 
here these models in dimensional form. The following terms have to be 
added to the source terms of Boussinesq equations when one wants to include 
bottom friction modelling. 

• Chezy law: 

= -^^^/^' 
where Cf is the Chezy coefficient. 
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Darcy-Weisbach law: 

Au lul 



' 8{h + 7]y 

where A is the so-called resistance value. This parameter is determined 
according to the simplified form of the Colebrook- White relation: 

^ = -2.03 log 



v/A ■ ^\UM{h + 7]) 

Here kg denotes the friction parameter, which depends on the compo- 
sition of the bottom. Typically kg can vary from 1mm for concrete to 
300mm for bottom with dense vegetation. 

Manning- Strickler law: 



2 U |U| 



4 1 
3 



{h + Tl) 

where k is the Manning roughness coefficient. 



5 Spectral Fourier method 

In this study we adopted a well-known and widely used spectral Fourier 
method. The main idea consists in discretizing the spatial derivatives using 
Fourier transforms. The effectiveness of this method is explained by two main 
reasons. First, the differentiation operation in Fourier transform space is ex- 
tremely simple due to the following property of Fourier transforms: /' = ikf . 
Secondly, there are very powerful tools for the fast and accurate computation 
of discrete Fourier transforms (DFT). So, spatial derivatives are computed 
with the following algorithm: 

1: /<-fft_(/) 

2: V ^ ikf 

3: /' ^ ifft iv) 
where k is the wavenumber. 

This approach, which is extremely efficient, has the drawbacks of almost 
all spectral methods. The first drawback consists in imposing periodic bound- 
ary conditions since we use DFT. The second drawback is that we can only 
handle simple geometries, namely, Cartesian products of ID intervals. For 
the purpose of academic research, this type of method is appropriate. 

Let us now consider the discretization of the dissipative Boussinesq equa- 
tions. We show in detail how the discretization is performed on equations 
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(22), (24). The other systems are discretized in the same way. We chose 
equations (22), (24) in order to avoid cumbersome expressions and make the 
description as clear as possible. 

Let us apply the Fourier transform to both sides of equations (22), (24): 



r/j = -ik ■ {h + er])u - -ht - —h^V^ht - ^hVh ■ Vht + —h^V^V ■ u 

e 2e e 6 

,,2 

+ ^K^V^hV -u + fx^h \Vhf V ■ u + fx^h^Vh ■ V(V ■ u), (50) 



U( + -£:?k|u|^ + ikrj + |k|^ u — -fi'^ikh'^V ■ Uj = 0, (51) 

where k = [kx, ky) denotes the Fourier transform parameters. 

Equations (50) and (51) constitute a system of ordinary differential equa- 
tions to be integrated numerically. In the present study we use the classical 
explicit fourth-order Runge-Kutta method. 

Remark on stability: A lot of researchers who integrated numerically 
the KdV equation noticed that the stability criterion has the form 

where A is the Courant-Friedrichs-Lewy (CFL) number and the number 
of points of discretization. In order to increase the time integration step 
At they solved exactly the linear part of the partial differential equation 
since the linear term is the one involving high frequencies and constraining 
the stability. This method, which is usually called the integrating factor 
method, allows an increase of the CFL number up to a factor ten, but it 
cannot fix the dependence on l/N"^. 

We do not have this difficulty because we use regularized dispersive terms. 
The regularization effect can be seen from equation (51). The same idea was 
exploited by Bona et al. [1981], who used the modified KdV equation (9). 

Let us briefly explain how we treat the non-linear terms. Since the time 
integration scheme is explicit, one can easily handle nonlinearities. For ex- 
ample the term [h + eri)u is computed as follows: 



{h + er])u = fft {{h + £ Re ifft (r/)) ■ Re ifft (u)) 



The other nonlinear terms are computed in the same way. 
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5.1 Validation of the numerical method 

One way to validate a numerical scheme is to compare the numerical re- 
sults with analytical solutions. Unfortunately, the authors did not succeed 
in deriving analytical solutions to the (1 + 1)D dissipative Boussinesq equa- 
tions over a flat bottom. But for validation purposes, one can neglect the 
viscous term. With this simpliflcation several solitary wave solutions can be 
obtained. We follow closely the work of Chen [1998]. In (1 -|- 1)D in the pres- 
ence of a flat bottom, the Boussinesq system without dissipation becomes 

r]t + u^ + e{ur])x - —u^xx = 0, (52) 

Ut + 'nx + £uux - -^^u^xt = 0. (53) 
We look for solitary-wave solutions travelling to the left in the form 

ri{x, t) = ?7(0 = vi^o + x + ct), u{x, t) = Br]{^), 

where we introduced the new variable ^ = Xq + x + ct and B, c, xq are 
constants. From the physical point of view this change of variables is nothing 
else than Galilean transformation. In other words we choose a new frame of 
reference which moves with the same celerity as the solitary wave. Since c is 
constant (there is no acceleration), the observer moving with the wave will 
see a steady picture. 

In the following primes denote derivation with respect to ^. Substituting 
this special form into the governing equations (52)-(53) gives 

ct]' + u' + eiurj)' — —u" = 0, 
6 

cu' + 7]' + euu' - c^u'" = 0. 
One can decrease the order of derivatives by integrating once: 

CT] + u + euT] u = U, 

6 

cu + ri -\ — u — c — u = U. 
' 2 2 

The solution is integrable on M and there are no integration constants, since 
a priori the solution behaviour at inflnity is known: the solitary wave is 
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exponentially small at large distances from the crest. Mathematically it can 
be expressed as 

lim r]{x,t) = lim u{x,t) = 0. 

^— >±oo 5^±oo 

Now we use the relation u{^) = -B?7(C) to eliminate the variable u from the 
system: 

(c + B)r]-B^7]" = -eB7]\ (54) 
6 

(1 + cB)r] - cB^r]" = -^B^V^. (55) 

In order to have non-trivial solutions both equations must be compatible. 
Compatibility conditions are obtained by comparing the coefficients of cor- 
responding terms in equations (54)-(55): 

Ib'--Bc = 1, 
2 2 

-B^-Bc = 0. 
6 

These relations can be thought as a system of linear equations with respect 
to B^ and Be. The unique solution of those equations is 

5 6 
Choosing i? > so that c > leads to 



'15 V15 

These constants determine the amplitude and the propagation speed of the 
solitary wave. In order to find the shape of the wave, one differentiates once 
equation (55): 

7r/'-/iV" = -IS^W- (56) 
The solution to this equation is well-known (see for example Chen [1998], 
Newell [1977]): 

Lemma 1. Let a, (3 he real constants; the equation 

has a solitary-wave solution if aP > 0. Moreover, the solitary-wave solution 
is 

r/(0=3«sech2fi./f(e + eo) 



where is an arbitrary constant. 
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Error between exact and numerical solutions 




100 150 200 250 300 350 400 450 500 
N - degrees of freedom 



Figure 8: Error on the numerical computation of a solitary wave solution. 
Here T = 1. 

Applying this lemma to equation (56) yields the following solution: 

r]{x,t) = -^sech^ I ^{x + ct + xq)] , (57) 



/ , 7VTE ,2/^7/ 
u[x,t) = — ^-^sech I "^(^ + + ^o) 

Note that this exact solitary wave solution is not physical. Indeed the 
velocity is negative whereas one expects it to be positive for a depression 
wave propagating to the left. In any case, the goal here is to validate the 
numerical computations by comparing with an exact solution. The method- 
ology is simple. We choose a solitary wave as initial condition and let it 
propagate during a certain time T with the spectral method. At the end of 
the computations one computes the L^o norm of the difference between the 
analytical solution (57) and the numerical one rj{x,T): 

€n ■■= max \r]{xi,T) - fj{xi,T)\ , 

l<i<N 

where {xi}^^^-^]^ are the discretization points. 
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parameter 


hi 


hr 


Xo 


Ax 


e /i z/i, z/2 


value 


0.5 


1.0 


-0.5 


0.3 


0.005 0.06 0.14 



Table 2: Typical values of the parameters used in the numerical computations 



Figure 8 shows the graph of as a function of A^. This result shows 
an excellent performance of this spectral method with an exponential con- 
vergence rate. In general, the error ejy is bounded below by the maximum 
between the error due to the time integration algorithm and floating point 
arithmetic precision. 

The exponential convergence rate to the exact solution is one of the fea- 
tures of spectral methods. It explains the success of these methods in several 
domains such as direct numerical simulation (DNS) of turbulence. One of the 
main drawbacks of spectral methods consists in the difficulties in handling 
complex geometries and various types of boundary conditions. 

6 Numerical results 

In this section we perform comparisons between the two dissipation mod- 
els (23) and (24). Even though the computations we show deal with a ID 
wave propagating in the negative x— direction, they have been performed 
with the 2D version of the code. The bathymetry z = —h{x,y) is chosen 
to be a regularized step function which is translated in the y— direction. A 
typical function h{x, y) is given by 

{hu X < Xo, 

hi + \{K. - hi) (l + sin {^{x - Xq - |Ax))) , Xq < x < Xq + Ax, 
hr, X > Xo + Ax. 

(58) 

This test case is interesting from a practical point of view since it clearly 
illustrates the phenomena of long wave reflection by bottom topography. The 
parameters used in this computation are given in Table 2. All values are given 
in nondimensional form. 

6.1 Construction of the initial condition 

We propagate on the free surface a so-called approximate soliton. Its clas- 
sical construction is as follows. We begin with the non-dissipative Boussinesq 
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equations on a flat bottom: 

r]t + ((1 + £v)u)^ - —u^xx = 0, 

Ut + Vx + -^{u )x - y^txxt = 0, (59) 
and look for u in the following form: 

u = -ri + eP + fj.^Q + 0{e^ + 6fi'^ + fj.^). (60) 

It is precisely at this step that one makes an approximation. One substitutes 
this asymptotic expansion into the governing equations and retains only the 
terms of order 0(e + /x^): 



m2 

r]t-Vx + ePx + l^'^Qx - 2er]r].^ + —Vxxx = 0{e^ + e^? + /), (61) 

6 



,2 



-Tit + + ePt + iJ^Qt + er/r/x + y^xxt = 0{e^ + efi^ + /). 
Add these two equations and set the coefficients of e and //^ equal to 0: 

e : Px + Pt- VVx = 0, (62) 

/i^ : Qx + Qt + IVxxx + IVxxt = 0. (63) 

Since the water depth is h = 1 + er] = 1 + 0(e), the approximate solitary 
wave should travel to the left with a celerity c = 1 + 0(e) and depend on the 
variable x + ct = x + 1 + 0{e) . Consequently one has the following relations: 

Replacing time derivatives by spatial ones in (62)-(63) yields 

Px "^Wxj Qx "^^xxx- 

By integration (using the fact that solitary waves tend to zero at infinity), 
one obtains 

P = lv^ Q = -lvxx (64) 
and the relation (60) connecting rj and u becomes 

u = -r] + ^T]^ - ^r]^x + .... (65) 
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Figure 9: Interaction between a left-running solitary wave and a step at two 
different times. The plots represent r]/10 while the true free-surface profiles 
are given hy z = erj. 



Substituting this expression for u into (61) yields a classical KdV equation 
for 77: 

3 \ 



- (^1 + ^sr]j r]^ - yVxxx = 0, (66) 
which admits solitary wave solutions of the form rj = r]{x + ct): 

r]{x,t) = ^^ sech^ (^-^^/6{c - l){x + ct)^ , 

where c > 1. The velocity u is obtained from (65) by simple substitution. 
This approximate soliton is used in the numerical computations. 



6.2 Comparison between the dissipative models 

The snapshot of the function r]{x,y,tQ) (divided by 10 for clarity's sake) 
during and just after reflection by the step is given on Figure 9. Recall 
that the free surface is given hj z = erj. Then we compare the two sets of 



6.2 Comparison between the dissipative models 



35 



Difference between model 2 and 1 at 1 = 0.4 

0.1 I < < < < < < 




X 



Figure 10: Free-surface snapshot before the interaction with the step: (left) 
the curves corresponding to the three models are almost superimposed; 
(right) difference between model II and model I. 

equations (22), (23) and (22), (24). To do so we look at the section of the 
free surface at y = along the propagation direction. 

Figure 10 shows that even at the beginning of the computations the two 
models give slightly different results. The amplitude of the pulse obtained 
with model I is smaller. It can be explained by the presence of the term uiSu 
which is bigger in magnitude than ei'2'V'^u. Within graphical accuracy, there 
is almost no difference between the conservative case and model II. 

In Figure 11 one can see that differences between the two solitons continue 
to grow. In particular we see an important drawback of the dissipation 
model I: just after the wave crest the free surface has some kind of residual 
deformation which is clearly non-physical. Our numerical experiments show 
that the amplitude of this residue depends almost linearly on the parameter 
1^1 . We could hardly predict this effect directly from the equations without 
numerical experiments. 

We would like to point out several soliton transformations in Figure 12 due 
to the interaction with bathymetry. First of all, since the depth decreases, the 
wave amplitude grows. Quantitatively speaking, the wave amplitude before 
the interaction is equal exactly to 8 (without including dissipation) and over 
the step it becomes roughly 9.4. On the other hand the soliton becomes less 
symmetric which is also expected. Because of periodic numerical boundary 
conditions we also observe the residue of the free-surface deformation coming 
through the left boundary. 

Figures 13, 14 and 15 show the process of wave reflection from the step at 
the bottom. The reflected wave clearly moves in the opposite direction. The 
fact that we see almost no difference between Model II and the conservative 
case should not lead to the interpretation that dissipative effects are not 



6.2 Comparison between the dissipative models 



36 



t= 1.100 




Figure 11: Free surface just before tlie interaction with the step 



t = 2.100 




■ No dissipation 



Model 1 

Model 2 



Figure 12: Beginning of the sohtary wave deformation under the change in 
bathymetry 
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t = 2.250 



No dissipation 




Figure 13: Initiation of the reflected wave separation 



t = 2.800 




Figure 14: Separation of the reflected wave 
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t = 3.200 

10 1 < < < < 




X 



Figure 15: Two separate waves moving in opposite directions 



important. One just has to wait long enough to see these effects play a role. 



7 Conclusions 

Comparisons have been made between two dissipation models. Model II, 
in which the decay is proportional to the second derivative of the velocity, 
appears to be better. At this stage we cannot show comparisons with lab- 
oratory experiments in order to demonstrate the performance of model II. 
Nevertheless, there is an indirect evidence. We refer one more time to the 
theoretical as well as experimental work of Bona et al. [1981]. In order to 
model wave trains, they added to the Korteweg-de Vries equation an ad-hoc 
dissipative term in the form of the Laplacian (but in ID). This term coincides 
with the results of our derivation if we model dissipation in the equations ac- 
cording to the second model. Their work shows excellent agreement between 
experiments and numerical solutions to dissipative KdV equation. Moreover 
our dissipative Boussinesq equations are in the same relationship with the 
classical Boussinesq equations [Peregrine, 1967] as Euler and Navier-Stokes 
equations. This is a second argument towards the physical pertinency of the 
results obtained with model II. 
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